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ABSTRACT 



We present a technique to synthesise telluric absorption and emission features both for in-situ wavelength calibration and for their 
removal from astronomical spectra. While the presented technique is applicable for a wide variety of optical and infrared spectra, we 
concentrate in this paper on selected high-resolution near-infrared spectra obtained with the CRIRES spectrograph to demonstrate its 
performance and limitation. We find that synthetic spectra reproduce telluric absorption features to about 2%, even close to saturated 
line cores. Thus, synthetic telluric spectra could be used to replace the observation of telluric standard stars, saving valuable observing 
time. This technique also provides a precise in-situ wavelength calibration, especially useful for high-resolution near-infrared spectra 
in the absence of other calibration sources. 
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1. Introduction 

Ground based astronomical observations suffer from the dis- 
turbing effects of the Earth's atmosphere. Strong ro-vibrational 
line systems of water vapour (H2O), carbon dioxide (CO2), and 
ozone (O3) shape the observational windows in the near in- 
frared that define the bandpasses of the J, H, K, L, and M fil- 
ters from 900 to 5500 nm. In addition, nitrous oxide (N2O) and 
methane (CH4) contribute noticeably to the atmospheric trans- 
mission losses in the near infrared, hampering the observation 
of important astrophysical lines. 

The Earth's atmosphere below about 60 km is in local ther- 
modynamic equilibrium and Kirchhoff's law applies. While the 
transitions of the molecules discussed above appear in absorp- 
tion in a stellar spectrum, they also appear in emission on the 
sky since the Earth's atmosphere is radiating its thermal en- 
ergy into space in these lines. The emission originating from 
the Earth's atmosphere competes with the astronomical signal 
at wavelengths longer than about 2300 nm, depending on ambi- 
ent temperature, slit width, and the readout noise of the detector. 
Because their spatial signature is extended and moderately fiat, 
beam-switching (nodding) techniques or in-situ fitting and sub- 
traction in 2D spectra is generally sufficient to remove the direct 
influence of the emission lines in the near and mid-infrared. 

The removal of the absorption features is considerably more 
difficult because no in-situ reference for these features is avail- 



* Based on observations made with ESO Telescopes at the Paranal 
Observatory under programme IDs 60.A-9076(A), 079.D-08 10(A), and 
179.C-015KA) 



able in the optical and near-infrared. The usual technique is to 
observe a star with either very few or well determined intrin- 
sic features (a "telluric standard") close in time and airmass 
to the scientific target. The choices for such a standard are 
early-type stars (spectral types mid B to late A) that exhibit a 
rather clean continuum despite strong hydrogen lines, or solar- 
like stars of spectral type early to mid G. I n the latter case , 
high resolution Fourier transf orm spectra (e.g jHase et al J, 12009: 
[Livingston & WallaceL 1 199 lh of the sun provide the necessary 
template to correct for the intrinsic stellar lines of the standard 
star. 

Despite the w ide use of this empirical calibration (e.g., 
IVacca et ail |2003[ and references therein), the use of standard 
stars for removing telluric absorption lines has several disad- 
vantages, especially when the correction is needed for low res- 
olution spectra. At a spectral resolving powers exceeding R ~ 
20000, telluric lines start becoming resolved and bandpasses 
with high telluric line density become accessible to ground based 
observations. The full potential of high resolution spectroscopy 
can only be used when unsaturated telluric lines are removed 
without compromising the intrinsic signal-to-noise ratio (SNR) 
of the data. To reach noise limited performance, the path length 
(airmass) and instrumental profile between the standard star and 
the science target must match. Also, the standard must be bright 
enough not to limit the final SNR. 

The intrinsic lines of early type standard stars, even when 
being accounted for, make the precise study o f near-infrared hy- 
drogen lines very difficult dVacca et al.L 120031) and the choice of 
a fast rotator can be an advantage. Also, the mismatches in the 
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Table 1. Log of CRIRES observations retrieved from the ESO archive. 



Object 


Spectral- 


Date & Time 


Atef 


V 


Slitwidth 


TTK 


Air- 


SNR C 




type" 




(nm) 


(GHz) 


(") 




mass 




HD 74575 


B1.5III 


2007-05-10 01:00 UT 


1134.0 


264,367 


0.2 


64,000 


1.30 


250 


HD 148605 


B3V 


2007-03-03 09:30 UT 


1504.8 


199,224 


0.2 


80,000 


1.03 


200 


HD 71155 


AOV 


2007-02-28 04:30 UT 


2076.5 


144,374 


0.4 


60,000 


1.24 


300 


HD 121847 


B8V 


2008-02-29 08:30 UT 


2336.2 


128,325 


0.4 


65,000 


1.01 


300 


HD 118716 


B1III 


2007-03-05 05:00 UT 


3907.8 


76,716 


0.2 


65,000 


1.30 


380 


HD 177756 


B9Vn 


2007-04-22 10:00 UT 


4760.8 


62,971 


0.2 


87,000 


1.06 


200 



a SIMBAD spectral type 

b Estimated resolving power from fitting telluric lines in the analysed spectral region. 
c Estimated signal-to-noise ratio in the continuum of the analysed spectral region. 



line depths of solar-type standard stars with the solar FTS atlas, 
caused, e.g., by differences in metallicity or effective tempera- 
ture, hamper precise abundance studies from near-infrared metal 
lines. Another problem arises for the L and M band, where the 
source flux of early- and solar-type standard stars can be the lim- 
iting factor, favoring the observation of non-stellar sources, like 
bright asteroids. 

Additional limitations in the use of empirical telluric cor- 
rections may be based on the spectrograph itself. Limited spec- 
tral coverage around the hydrogen lines of early-type standard 
stars prohibit fitting and correcting these lines, which leaves 
a disturbed continuum in the spectrum of the science target. 
Adaptive optics fed long-slit spectrographs, such as CRIRES at 
the VLT (see Sect. [2]), introduce another level of complication. 
The diffraction-limited full width at half maximum (FWHM) 
of the stellar point spread function (PSF) can be considerably 
smaller than the slit width, so that one often observes with an 
under-filled entrance slit. The actual instrumental profile, and 
hence spectral resolving poweiQ, is then governed by the PSF 
of the source, rather than by the slit width of the spectrograph. 
Moreover, slight displacements of the source in the slit will lead 
to shifts in the wavelength at the detector plane. In the case of 
CRIRES, a shift of only 30 m/s (1/50 of a CRIRES pixel or 4 mas 
at the entrance slit) can notably degrade the telluric correction 
(see, e.g., Sec. 14. 1.2b . Since most of the telluric lines in the near- 
infrared are not fully resolved, even at /?=100000, their shapes 
are governed by the instrumental profile. As a consequence, even 
at the same airmass the telluric standard star cannot perfectly 
correct the telluric features in the spectrum of the science tar- 
get because differences in the performance of the adaptive op- 
tics ultimately lead to differences in the instrumental profile and 
variations of the wavelength solution. A further drawback in the 
use of telluric standard stars is the waste of observing time for 
the acquisition and observation of these stars, especially under 
the constraint of matching airmasses at small time differences. 
As mentioned earlier the problem of source brightness becomes 
especially apparent at longer wavelengths where one can easily 
spend as much or even more time observing a telluric standard 
than a science target. 

While telluric emission and absorption features have been 
discussed so far as a disturbance in astronomical spectroscopy, 
they can as well turn out to be very useful. Determining the 
wavelength solution of near infrared spectra, especially at high 
spatial and spectral resolution, proves to be very difficult. Rare 
gas emission line lamps used in low resolution near-infrared 
spectrographs, such as He, Ne, Xe, and Kr lamps, provide much 

1 Throughout this paper we define the spectral resolving power as 
R = A/AA, where AA is the FWHM of the instrumental profile. 



too sparse of a line coverage to be useful at high spectral reso- 
lution. Even the rich spectra of ThAr emission line lamps, com- 
monly employed by high resolution optical spectrographs, have 
a much lower line density in the near-infrared than in the red 
optical. Typical line densities of ~ 400 lines per AA = 100 nm 
around A = 1000 nm drop quickly to less than ~ 20 l ines per 
AA = 100 nm at A = 2500 nm dKerber et all 120081) . Given 
the typical wavelength coverage of ~ A/200 per detector, many 
CRIRES spectral settings remain poorly calibrated. Moreover 
the dynamic range of the ThAr emission spectrum is very high, 
owing to the strong contrast between the Th and Ar lines, leav- 
ing weaker lines at low SNR while nearby stronger lines saturate 
quickly. 

Recently frequency combs (lUdem et all 120021) have been 
proposed for the ideal spectral calibration of high resolution 
spectrographs. In principle a frequency comb would provide for 
high line density while linking the dispersion solution of the 
spectrograph to the time standard. However, the line density of 
a simple comb is much too high, so that some optical filtering 
is required, which will produce detrimental artifacts. Also, a fre- 
quency comb might solve for the problem of insufficient line 
density and uniformity of emission line lamps but the instru- 
mental profile will still be different than for the science target, 
ultimately limiting the final precision of the wavelength calibra- 
tion of the science spectrum. 

This limitation applies also for atmospheric emission lines 
from ro-vibrational levels of the Hydroxyl radical OH* (Meinel 
bands, commonly referred to as OH* airglow emission) which 
are routinely used for the c alibration of low-reso lution near- 
infrared spectrographs (e.g., Rousselot et all l2000h . Moreover, 
neither the line strength nor the line density is sufficient to be 
used as a useful wavelength calibration over a wide spectral 
range at high spectral and spatial resolution, as provided by 
CRIRES. The small spatial pixel scale of CRIRES (86mas/pix) 
makes the use of OH* lines for calibration purposes challenging, 
given that each pixel is less than 10~ 2 arcsec 2 on the sky. In con- 
trast, the atmospheric absorption and emission features of other 
molecules show a high density of strong absorption lines and 
thus provide a natural in-situ wavelength calibration. Moreover, 
the telluric absorption lines suffer the same instrumental effects 
as the lines in the stellar spectra and, hence, provide an intrin- 
sically more precise calibration source as compared to telluric 
emission lines or lamp emission lines, the latter even recorded in 
separate exposures and often hours apart from the science spec- 
tra. 

We have thus investigated the feasibility of calculating theo- 
retical atmospheric transmission spectra in order to use them for 
wavelength calibration and to remove their signature from the as- 
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tronomical spectra, thus, replacing standard stars as a source of 
telluric calibrato rs. This idea is not new and has been used suc- 
cessfully before (|Lallement et al.L Il993t IWidemann et all 1 19941 
Baile y et aD.l2007l) . Nevertheless, modeling telluric spectra has 
not evolved to a standard technique in optical and infrared spec- 
tral data reduction because suitable radiative transfer codes were 
hard to access and synthesising adequate spectra was hindered 
by the incompleteness of molecular line databases. This situa- 
tion has improved in the last years. 

In an effort to characterise the abilities and the efficient us- 
age of the CRIRES spectrograph, we have developed a general 
method to model the transmission and emission spectrum of the 
Earth's atmosphere above Cerro Paranal. This method can be 
used to perform a wavelength calibration of CRIRES and to sub- 
sequently remove the telluric absorption lines apparent in these 
spectra. Moreover the model can also be used to predict the per- 
formance of planned CRIRES observations using a forecast of 
the atmospheric conditions during a planned observing run. 

The paper is structured as follows: In Sect. [2] we describe the 
data reduction of the CRIRES observations. In Sect. [3] we de- 
scribe the the radiative transfer code and input parameters used 
to synthesise the telluric spectra, which we compare to the ob- 
servations in Sect. [4] We summarise our findings in Sect. [5] and 
close with conclusions about the performance of the presented 
technique. 



2. Observations and data reduction 

CRIRES is a high resolution (R < 100000) near-infrared 
(/Ll=960-5200nm) adaptive opti cs (AO) - f ed spectrograph 
at the VLT on Paranal, Chile (iKaufl et all 12004 |2006a b; 
iPaufique et all l2006h . Commissioning and science verification 
(SV) observations were executed in October 2006 and February 
2007. The spectrograph is available to the community and has 
been used for regular observations since April 2007. 

We selected observations from both SV periods as well as 
more recent observations from the ESO data archive to test the 
telluric model spectra described in Sect. [3] For this purpose we 
chose observations of standard stars at different wavelengths, 
observed at different airmasses with nominal spectral resolving 
powers of R = 50000 and 100000. The log of observations is 
given in Tab.Q] 

Observations were always obtained in an AB or ABBA nod- 
ding pattern. Data reduction followed the standard steps for 
long-slit spectrographs as applicable to the CRIRES data for- 
mat. All raw frames were treated with a non-linearity correction 
before pairwise subtraction removed the atmospheric emission 
features. The individual A-B and B-A frames were then divided 
by a normalised flatfield. Optimally extracted ID spectra at both 
nodding positi ons were obta ined by using a custom made IDL 
script based on lHornd d 1986b . 

Due to the curvature of the slit, it is inadvisable to combine 
the 2D spectra before the extraction and wavelength calibration 
step as is done in the current version of the CRIRES data re- 
duction pipeline. The amount of curvature can be of the order of 
1 pixel between two nodding positions at a nod throw of Wand 
is variable in both spatial and spectral directions. We thus com- 
bined the individual ID spectra after correcting for the slit cur- 
vature by producing a telluric model spectrum using the method 
outlined in Sect. [3] and interpolating the wavelength calibrated 
spectra to a common wavelength vector. All reduction steps were 
performed individually for the four chips of the CRIRES focal 
plane array. 



3. Spectral synthesis 

To compute theoretical transmission and radiance spectra of the 
Earth's atmosphere, a radiative transfer code is used that takes as 
an input a model of the atmosphere in terms of vertical temper- 
ature, pressure and molecular abundance profiles, as well as line 
data (frequency, line strength, pressure broadening coefficients, 
etc.). In this section we describe the usage of such a code as well 
as the necessary inputs in further detail. 

3.1. Radiative transfer code 

There are several radiative transfer codes tailored for the con- 
struction of telluric spectra. These codes are based on line-by- 
line computations o f a la y ered m odel of the Earth's atmosphere. 

STRANSAC (IScottL 1 1974 is one the oldest of such 
codes. It utilises the GEISA spectral line catalogue 
( Jacqu inet-Husson et al l [T999[l200l . 

The 4A (Automatized Atmospheric Absorption Atlas) code 
also uses GEISA, and it provides a fast and accurate line-by-line 
radiative transfer model that is particularly efficient in the in- 
frared region between 3,000 nm and 17,000 nm. The latest ver- 
sion of 4A can be retrieved free of charge for non-commercial 
usage from NOVELTISfl France. 

ATRAN dLordl 1 19921) is another code to compute synthetic 
spectra of atmospheric transmission. It is based on a fixed atmo- 
spheric layering (see Sect. 13.21 ) and uses the HITRAN database 
(see Sect. 13 .3b for line data input. 

RFM (Reference Forward Model) is a line-by-line radiative 
transfer model originally developed at Oxford UniversitjQ, under 
an ES A contract to provide reference spectral calculations for the 
MIPAS instrument launched on the ENVISAT satellite in 2002. 

The code uses HITRAN 2000 as its line database. 

In t his pa per we only consider FAS CODE (Clou gh et all 

fl98ll fl992b . and used both its commercial version 
FASCODE3P/PcLnWirQ and the free of charge version 
LBLRTMB LBLRTM is available as FORTRAN source code 
and can be compiled on various platforms. Its most common ap- 
plication is the spectral retrieval from ground and satellite based 
measurements in climatology. The code accepts a custom model 
atmosphere and uses the HITRAN database (see Sect. 13.31 ) as 
input for the line data. We found LBLRTM easy-to-use, yet 
flexible, and computed all spectra presented in this paper with 
this code. 

3.2. Input I: Model atmosphere 

The choice of a representative model atmosphere is crucial for 
a realistic theoretical spectrum of the atmosphere. One of the 
more basic descriptions of a model atmosphere is given by the 
1976 US Standard Atmospher^ and its geographical subdivi- 
sions (tropical, mid-latitude and sub-arctic summer and win- 
ter models). Due to the geographical location of Cerro Paranal 
and the fractional changes in the abundance of trace gases since 
its definition, the 1976 US Standard Atmosphere gives a rather 
poor representation of the atmospheric conditions above Cerro 
Paranal. For example, the average surface temperature of the 

2 http://www.noveltis.net/4A0P/ 

3 http://www.atm.ox.ac.uk/RFM/ 

4 http : //www . ontar . com/Sof tware/Products . aspx 

5 http : //rtweb . aer . com/lblrtm_description. html 

6 http : //ntrs . nasa . gov/archive/nasa/casi . ntrs . nasa . gov/ 
19770OO9539.1977OO9539.pdf 
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Fig. 1. Vertical profiles of the 10 most abundant molecules in- 
cluded in our model atmosphere. The water vapour (H2O) abun- 
dance for the lowest atmospheric layers (surface height < 26 km) 
is replaced by a meteorological model (GDAS or MM5) for the 
time of observation. 



model at 2640 m is only 271 K and thus too low by about 10 K 
on average. 

In addition to the semi-static distribution of atmospheric 
trace gases, one has to account for the slowly variable temper- 
ature and pressure levels as well as the highly variable water 
vapour content. Thus, one has to base the calculations on a more 
flexible and realistic model atmosphere which takes meteorolog- 
ical inputs into account. 

We consider here the use of meteorological models from the 
Air Resources Laboratory (ARL) at the US National Oceanic 
and Atmospheric Administration (NOAA). Sounding files for 
temperature, pressure and dew point temperature in the tropo- 
sphere and lower stratosphere (surface height < 26 km) are avail- 
able at the ARL websiteQbased on either Global Forecast System 
(GFS) or archived Global Data Assimilation System (GDAS) 
models. GFS analysis data are available in 3 hr intervals for the 
next 84 hours and 12 hr intervals for the next 384 hours. GDAS 
models are available in 3 hr intervals for all dates since Dec 
2004. 

The horizontal resolution of the GFS and GDAS models is 
1° (approx. 110km). Given that Cerro Paranal is very close to 
the Pacific ocean (~12km), the influence of the ocean climate 
on the chosen grid point could compromise its validity. We have 
thus also computed model atmospheres for three nights during 
CRIRES S V periods on October 10, 2006, and March 3-4, 2007, 
based on the Fi fth-Generat i on NC AR / Penn State Mesoscale 
Model (MM5) dGrell et all 1 19951) . which uses the operational 
analyses of the European Centre for Medium-Range Weather 
Forecasts (ECMWF) as start and boundary conditions and was 
refined for the local topography of the Paranal area with a final 
resolution of 1 km in the inner integration area. Both meteoro- 
logical models reproduced the ambient pressure and temperature 
at the VLT as recorded in the file headers of the analysed obser- 
vations to within 5 hPa and 1 .5 K, respectively. 

All atmospheric models have to be supplemented with infor- 
mation on the vertical distribution of all molecules other than 
H2O and used in the radiative transfer calculations. Moreover, 
since the meteorological models are limited to estimated surface 
heights of ~26km, we also have to add temperature, pressure 



and water vapour profiles for the remaining atmospheric layers 
(26 < h < 75 km). For this purpose we have chosen an equato- 
rial MIPAS model atmospherqj, constructed by John Remedios 
(U. Leicester). This model contains temperature, pressure and 
abundance information for 30 molecular species in 121 levels 
from 0-120 km at 1 km spacing. In Fig. Q]we show the mixing 
ratios of the 10 most abundant molecules in this model. The ac- 
tual amount of molecules per cm 3 drops roughly exponentially 
with the surface height and the predominant part of the telluric 
transmission is shaped for most species in the first few km above 
the observatory. 

In Fig. [2] we compare our modified MM5 model to a GDAS 
model for March 3, 2007 at 09:00 UT. The MM5 was u sed with 
a set of modifications developed by Zangl (2002, 120031) in order 
to improve the numerical accuracy over steep topography. While 
temperature and pressure profiles are nearly identical, the water 
vapour profile is noticeably different. The GDAS model predicts 
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Fig. 2. Comparison of atmospheric models for Cerro Paranal, 
GDAS (red) and MM5 (black). Shown are pressure (upper left 
panel), temperature (upper right panel) and water vapour (lower 
panel) profiles for March 3, 2007 at 09:00 UT for surface heights 
above Paranal (2.64km). Information for heights greater than 
26 km (GDAS) and 22 km (MM5) were delivered by a global 
model for equatorial regions. 

a total water vapour column, PWV0, of 4. 1 mm, the modified 
MM5 model predicts 2.4 mm. Moreover, the vertical distribu- 
tion of the water vapour is different as well. The GDAS model 



http : //www . arl . noaa . gov/ready/cmet . html 



8 http : //www-atm . physics . ox . ac . uk/RFM/ atm/ 

9 measured as the precipitable water vapour amount at zenith 
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predicts a much more humid surface layer, but dries out beyond 
the prediction of the MM5 model above 6 km height. A similar 
situation was found for October 10, 2006. Both models predict 
the same PWV of 1 .8 mm but the GDAS model is again more hu- 
mid below 3 km and less humid above. In any case, the models 
are often too wet when compared to the observations. The total 
water vapour amount was always scaled by 40%-70% to match 
the observed line depth (see §4). The different vertical distribu- 
tions of the water vapor did not lead to a significant difference 
in the fit quality of the water vapor lines once the total water va- 
por amount was rescaled to reach the best fit. We compare the 
predicted and fitted PWV values in Tab. [2] and conclude here, 
that the GDAS models are applicable and sufficient as meteoro- 
logical input to refine the global equatorial model atmosphere. 
See Sect. [4] for a further discussion about the performance and 
limitation of the actual fitting. 

3.3. Input II: Molecular line database 

HITRAN, in its latest 2008 edition (iRothman et all |2009), is a 
molecular line database containing frequency, line strength and 
pressure broadening coefficients for more than 1 .7 million spec- 
tral lines of 42 different molecules and their common isotopes. 
The database is regularly updated and supplemented with new 
line data and molecular cross-sections. LBLRTM is distributed 
with a line file that is based on the HITRAN 2004 version and in- 
cludes individual updates for some molecules until June 2009. A 
complete update to HITRAN 2008 is planned for the near future 
(Mark Shephard, private communication). 

4. Results 

We have computed synthetic telluric transmission spectra for 
a number of selected CRIRES observations (see Tab. [TJ using 
LBLRTM and a GDAS model atmosphere closest in time to the 
observation. The wavelength range of the model was set to the 
coverage of the chosen CRIRES setting. LBLRTM uses the al- 
titude angle of the telescope pointing, the height of the obser- 
vatory and the model atmosphere to calculate the line-by-line 
and layer-by-layer transmission and radiance. We make no use 
of the calculated radiance at this point, since the telluric emis- 
sion is removed conveniently by nodded observations, and use 
the transmission spectrum only. 

An IDL based Marquard t-type x 1 minimisation code 
dBevington & Robinso r3. [T99l p. 161) was used to (a) refine 
the wavelength solution by fitting a second order polynomial, 
using the original wavelength calibration as first guess; (b) re- 
determine the total water vapour amount and column density of 
other molecules present in the chosen setting; and (c) determine 
the instrumental profile and continuum rectification parameters. 

The determination of the instrumental profile (IP) is crucial 
for a proper fitting of the only partially resolved atmospheric 
lines. We have used a singular- value de composi t ion (S VP) tech- 
nique in an implementation given by iRucinskil d 1999b to deter- 
mine the instrumental profile. The wavelength scale was interpo- 
lated to a five times higher sampling than the original, Nyquist 
sampled wavelength scale and a 20 pixel wide instrumental pro- 
file was constructed from 11-13 Eigenvalues from all telluric 
lines present in the spectrum. Higher number of Eigenvalues did 
not lead to a further improvement of the fit but amplified the 
noise in the instrumental profile. In all cases, the instrumental 
profile was extremely well matched by a single Gaussian func- 
tion whose FWHM depends on the the slitwidth and the perfor- 
mance of the adaptive optics, see Fig. [3] 



In cases where extremely narrow spectral lines were present 
(e.g. the unresolved ozone lines at 4700 nm - see Sect. 14. 1 .51 up 
to 30 Eigenvalues could be used without amplifying the noise. 
This analysis confirmed that a Gaussian model is a valid approx- 
imation for the instrumental profile at signal-to-noise levels be- 
low ~ 300. We have thus used a single Gaussian model for the 
final fitting of the telluric model, the width of the profile being 
the only free parameter. We note that our determined FWHM are 
broader than expected for R = 100 000. At the time most of the 
data presented in this paper were recorded, the focal plane array 
of CRIRES was slightly out of focus. Thus, the maximum re- 
solving power measured in the presented dataset is R ~ 87 000. 

4. 1 . Performance in selected settings 

We present CRIRES observations of bright standard stars that 
contain clear and isolated line systems of the most abundant 
molecules of the Earth's atmosphere. We evaluate the perfor- 
mance of the telluric correction by checking the residuals for 
systematics by a rigid quantisation of the signal-to-noise ratio 
(SNR) in the corrected spectra. We chose the SNR as a mea- 
sure for the performance of the telluric model fit since it has 
the strongest practical implication for the observer, as it allows 
a judgement as to which precision a position or an equivalent 
width of a stellar line can be measured in a spectrum. In regions 
of telluric contamination, the residuals from the removal of tel- 
luric lines are the main factor hindering a high signal-to-noise 
measurement, otherwise limited only by the source flux or the 
uniformity of the detector response. 

In order to compare the SNR before and after the correc- 
tion, we first compute the intrinsic noise level of our CRIRES 
spectra by adding the photon noise of the source and the mea- 
sured background noise in quadrature for each pixel along the 
dispersion direction. The background noise consists of the de- 
tector readout noise, the photon noise of the dark current and the 
photon noise of the thermal emission of the atmosphere and the 
telescope after pairwise subtraction of nodded frames. We have 
validated this procedure by measuring the standard deviation in 
the continuum level of each spectrum. The result was always 
consistent with our noise estimate. 
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Fig. 3. CRIRES instrumental profile, obtained from SVD fit- 
ting of methane (CH4) lines in setting /l re f 2336 nm (see Fig.|6]l. 
Black: empirical profile constructed from 21 Eigenvalues. Red: 
best fit Gaussian model with a FWHM of 2.67 pixel (4.0km/s, 
fl=70,000). 
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Hence, for the input spectra we have the SNR as a function 
of wavelength which scales with the square root of the received 
flux, and thus with the square root of the atmospheric transmis- 
sion T. In addition, the telluric emission lines in the thermal in- 
frared example (A ~ 4700 nm, Sec. I4.1.51 l add noise where the 
lower transmission of the same telluric lines in absorption de- 
grade the SNR in comparison to the continuum level. 

The fitting procedure minimises the residuals in the 
observed-computed spectrum (O-C) over the whole wave- 
length region of the considered spectral range of ~ /i/200, 
recorded on one CRIRES detector. For the removal of the tel- 
luric lines, however, the measured spectrum has to be divided by 
the computed model. After division by the model, we measure 
the noise in the spectrum as the standard deviation of all pixels. 
The SNR is then simply the inverse of the measured noise, given 
that after the division the signal level is at unity for a normalized 
spectrum. 

It is important to remember that even for a perfect model of 
the telluric lines the O/C residuals will be larger, and the SNR 
smaller, in comparison to the undisturbed continuum, since the 
spectrograph receives less flux where telluric absorption lines 
are present. Hence, when dividing through the telluric transmis- 
sion spectrum, and thus correcting the received flux for the pres- 
ence of these lines, the noise in these regions is higher than in the 
neighboring continuum, simply because less flux was received in 
the lines in the first place. Albeit this caveat, we use the formal 
SNR as a direct measure of the performance of the modeling 
since it reflects the impact of the telluric lines on the astronom- 
ical spectrum and gives the observer a meaningful number to 
work with. 



4.1.1. H 2 at 1505nm 

Water vapour (H2O) is the most important absorber in the near- 
infrared. In addition to the strong v\ and V3 fundamental bands 
at 2700 nm, several overtone and combination bands centered at 
940 nm, 1100 nm, 1380 nm, 1870nm, and 3200 nm are present in 
near-infrared spectra. We picked a few isolated, unsaturated H2O 
lines at 1505 nm in a CRIRES standard star setting from March 
03, 2007. The airmass of the observation was 1.033, the resolv- 
ing power was -80,000 at a slitwidth of 0.2". Both the GDAS 
model and the modified MM5 model had to be scaled to a PWV 
of 1.6 mm for the best fit with the observation (see Fig. |4j». The 
achieved correction using the GDAS and MM5 models is gen- 
erally satisfying but falls behind the expectation from the noise 
limit. The average SNR of the measured spectrum is ~200. After 
the division by the model spectrum, the average SNR drops to 
~180. In the line centers (40-70% transmission) the SNR drops 
from ~160 to ~50 after dividing by the model. Note how lines of 
equal strength, e.g., around 1504 nm are fitted with greatly vary- 
ing precision, pointing to uncertainties in the line parameters. 

4.1.2. C0 2 at 2060 nm 

Carbon dioxide (CO2) has a nearly constant vertical mixing ra- 
tio (see Fig. [1} but due to its strong interaction with the bio- 
sphere, the low altitude CO2 abundance can be expected to be 
variable in time. CO2 is the third most abundant molecule after 
oxygen and nitrogen at heights exceeding ~4km and exhibits 
various strong absorption bands c entered at 1400 nm, 1600nm , 
2000 nm, 2700 nm, and 4300 nm. iKenworfhv & Hanson! |2004) 
explored the performance of empirical telluric line removal 
around the important He I lines at 2058.69 nm, which is buried 




1504 1505 
Wavelength (nm) 



Fig. 4. Water vapour (H 2 0) lines on chip three of CRIRES set- 
ting /t le f 1504.8 nm. Top panel: CRIRES spectrum (black), and 
synthetic transmission spectrum (red). Lower panels: O-C and 
O/C residuals of the fit. 

under a comb of lines from the R branch of the 200 1 3 <— 00001 
transition of telluric 12 C 16 16 0. 

We have used a standard star observation obtained on March 
01, 2007 in this spectral region to test the performance of our 
synthetic telluric model. The airmass of the observation was 
1.24, the resolving power was ~60,000 at a slitwidth of 0.4". 
The CO2 abundance is given by the global atmospheric model 
(see Sect. 13.21 ) and had to be scaled down by ~1% to get the best 
match with the observation (see Fig. [5j, despite the fact that the 
atmospheric CO2 abundance is continuouslygrowing and shows 
yearly changes of only + 1 % (NOAA/ESRL0). 




2058 2060 2062 2064 2066 2068 

Wavelength (nm) 

Fig. 5. Same as Fig. [4] but for carbon dioxide (CO2) lines on 
chip two of CRIRES setting A re f 2076.5 nm. Regions with tel- 
luric transmission of less than 5% are not shown in the residuals. 



The average SNR of the spectrum is ~300. It drops to about 
230 at transmissions > 60% when corrected for the telluric ab- 
sorption. The SNR at the line cores (10-50% transmission) is 
~160 and drop s to ~50 after the correction, see Fig. [5] As al- 
ready noted by IKenworfhv & Hanson! d2004l) . a precise wave- 
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length solution is crucial for a good fitting. Wavelength errors at 
the level of 1/50 pixel (~30m/s) notably degrade the fitting qual- 
ity because of the steep line profiles. While the line positions in 
the HITRAN database have uncertainties at the m/s level, the un- 
certainty in the pressure shifts may become an issue at this level 
of precision. 

After dividing by the best fit synthetic model, the region 
around the He I line has a SNR of >160. We note that a high- 
/ line of 13 CC>2 leaves a 3% residual right at the rest wavelength 
of the He I line. We also note that for increasing J quantum num- 
bers, the line shape of the 12 CC>2 is less well matched as close to 
the band center. CO2 exhibits strong line coupling which is taken 
into account by LBLRTM (see discussion in Sect. [5]). It remains 
thus unclear whether the systematic increase of the residuals are 
due to uncertainties in HITRAN or due to insufficient treatment 
of the line coupling. 

The 12 CC>2 lines are fully resolved and only a few narrow 
lines are present in the setting. This inhibits a precise determina- 
tion of the instrumental profile, which becomes very sensitive to 
noise in the spectrum. Thus, the width of the instrumental profile 
is partly degenerate with the line strength of CO2 in the model. 
The exact locus of the best correction can be modified by slight 
changes of the FWHM of the instrumental profile. This favours a 
customised telluric correction, tailored for a specific stellar line. 

4.1.3. CH 4 at 2320 nm 

Methane (CH4) exhibits strong bands centered at 1660nm, 
2200 nm, and between 2300 nm to 3830 nm, although it is 
much less abundant than H2O and CO2. The methane abun- 
dance in our model is set by the global atmospheric model 
(see Sect. I3.2l >. Methane is another of the so called greenhouse 
gases. Its atmospheric concentration is slowly increasing with 
time, but its growth rate d ropped considerably in the last years 
dDlug okenckv et all I2003I) . We chose a observation from Feb 
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\ 1.00 
o 
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Fig. 6. Same as Fig. |4] but for methane (CH4) and a few 
water vapour (H2O) lines on chip two of CRIRES setting 
A Te f 2336.2 nm. 

29, 2008, covering a strong CH4 band at 2320 nm for fitting. The 
airmass of the observation was 1.01, the resolving power was 
~65,000 at a slitwidth of 0.4". This region is populated by CO 
first overtone lines in astrophysical objects and thus marks a very 
interesting spectral region. Besides the strong CH4 lines, there 
are also three prominent water vapour lines in the covered re- 



gion, one strong line at 2321.95 nm, and two weaker line blends 
at 2322.66 nm and 2323.27 nm. The best fit PWV at zenith was 
2.2 mm. We did not have to scale the original CH4 abundance in 
the equatorial model. The average SNR of the spectrum is ~300, 
dropping to -220 when corrected for the telluric absorption. 

4.1.4. N 2 Oat3907nm 

An observation from March 05, 2007 covered a nitrous oxide 
(N2O) band around 3907 nm at airmass 1 .3. The resolving power 
was only ~65,000 at a slitwidth of 0.2". N2O is a factor of two 
less abundant than CH4. Still, it exhibits strong line systems be- 
tween 2870 and 4500 nm. As for C0 2 , we had to scale the N 2 
abundance in our model down by ~5% to match observations. In 
addition to the strong N2O lines, there are weaker lines of CH4 
and H2O present in this setting. The best fit PWV at zenith was 
8.6 mm. The average SNR of the spectrum is -380, dropping 
only by 20% to ~3 10 when corrected for the telluric absorption. 




f\f) 



3905 3910 
Wavelength (nm) 



Fig. 7. Same as Fig.|4]but for nitrous oxide (N2O) lines on chip 
three of CRIRES setting A re f 3907.2 nm. Minor contributions be- 
tween AA 3898-3905 nm arise from H 2 Q and CH 4 . 



4.1 .5. H 2 0, 3 , C0 2 , and CO at 4700 nm 

Ozone (O3) is a major opacity source in the UV where the 
Hartley and Huggins bands block most of the stellar and solar ra- 
diation. The Chapp uis ban ds at 600 nm are weak, but still visible 
in optical spectra ( GriffinL [20051) and originate from the strato- 
sphere where ozone has its maximum abundance at about 30 km 
height (see Fig. [TJ. The near-infrared band of O3, centered at 
4750 nm, is produced by overtone and combination transitions. 
These lines are strong but comparably narrow, which is due to 
the low atmospheric pressure at the height of formation. 

We chose a standard star observation from April 22, 2007 
covering a region around 4700 nm at airmass 1 .06. The resolving 
power was ~87,000 at a slitwidth of 0.2". In the chosen region, 
H2O, CO2, and CO lines contribute to the opacity, see Fig. [8] 
The Q branch of the 20001 <-01 101 transition of C0 2 is very 
prominent (see the second to last panel in Fig. [8j. 

The fundamental transitions of carbon monoxide (CO) are of 
great importance for astronomical observations. Unfortunately, 
the CO lines are also very strong in the telluric spectrum and 
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the transitions with 2 < J < 12 are saturated, even at mod- 
est airmasses. High spectral resolution and velocity offsets from 
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Fig. 8. Transmission of the individual species involved in the 
spectral region shown in Fig. [9] From top to bottom: water 
vapour (H2O), ozone (O3), carbon dioxide (CO2), and carbone 
monoxide (CO). The transmissions are shown before convolu- 
tion with the instrumental profile. 



barycentric motion and radial velocity of the science target are 
needed to disentangle the stellar and telluric components. The 
average SNR of the measured spectrum in Fig|9]is ~200, which 
dropped to ~100 when divided by the model spectrum. The CO2 
Q branch is not well fitted and peaks out as a strong residual 
near 4695 nm. The second strongest residual near 4707.5 nm is 
a H2O line. This line was much better fitted with the line data 
from the HITRAN 2004 database before the update in 10/2008 
(AER V2.1). Its line strength is notably different in the newest 
line file (AER V2.2) and causes a worse fit. Other minor resid- 
uals are mostly from O3. Because the line formation happens in 
a cold part of the Earth's atmosphere at ~20-30km height and 
its abundance sharply decreases below 15 km, O3 does not con- 
tribute notably to the telluric emission received at the ground 
but is a strong source of absorption. Only small scaling correc- 
tions at the order of 2.5%-3% had to be applied to the model 
abundances of CO2, and CO while the O3 abundance had to be 
increased by 6%. The best fit PWV was 7.2 mm (7.7 mm using 
the HITRAN 2004 database). 

4.2. Limitations 

For the cases presented in the previous subsection, we have never 
reached the formal noise limit of the observations but we found 
a satisfying fit of the model to the data. To illustrate the limits 
of our approach, we present an extreme example of very badly 
fitting H2O lines in Fig. [10] The airmass of the observation was 

1.3, the resolving power was ~64,000 at a slitwidth of 0.2". This 
is the worst fit we have come across and it is clearly related with 
the line strength information in HITRAN. Neighboring lines in 
the model are both too strong and too weak. The predicted PWV 
at zenith was ~2 mm and was scaled down during the optimi- 
sation to 0.35 mm. We find this reasonable, since the PWV was 
reported at the same level for other spectral settings taken during 
the same night where the water vapour lines were fitted much 
better. The lines in the setting shown in Fig. [10] are thus very 
strong lines as they show up with an average transmission of 
50% and below even under such dry conditions. Water vapour 
lines also caused problems in the setting presented in Sec. 14. 1.51 
and left residuals in other settings. It is thus the water vapour 
which appears the hardest to model and to be affected the most 
from uncertainties in the line database. 



4695 4700 4705 4710 

Wavelength (nm) 



Fig. 9. Same as Fig.[4]but for ozone (O3), water vapour (H2O), 
carbon dioxide (CO2) and carbone monoxide (CO) on chip one 
of CRIRES setting ^ ref 4760.8 nm. 




1132 1133 1134 1135 1136 

Wavelength (nm) 



Fig. 10. Same as Fig.|4]but for water vapour (H2O) lines on chip 
three of CRIRES setting A ref 11 34.0 nm. 
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4.2.1. Limitations from the HITRAN database 

The HITRAN database contains estimators on the uncertainty 
of p osition, line strength and pressure broadening for all en- 
tries (Rothman et alll2005l) . For example, the 12 CC>2 lines shown 
in Fig. [5] have an average uncertaint y in their positions of ~ 
5 x 1(T 5 cm -1 (Miller & Brownl l2004T) or ~ 3 m/ s and the uncer- 
tainty in the line strength is at the few-% level ( Ta shkun et all 
2003). The strongest CH4 lines shown in Fig.|6]have an average 
uncertainty in their positions of ~ 10 3 cm" or ~ 60 m/s and un- 
certainties in their hne s trengths of up to 10% ( Roth man et all 
120051; iBrown et al l l2003l) . This level of positional uncertainties 
influences the wavelength solution and the line fitting. As ex- 
pected, even the strongest water vapour lines shown in Fig. [10] 
have positional uncertainties of up to 0.01 cirT 1 or ~ 2.6 km/s 
and their line s trengt h can be uncertain to more than 20% 
dRofhman et all 120051) . This explains the bad fit to the data 
shown in Fig.fTOl 



4.2.2. Limitations from LBLRTM 

Other limitations arise from the limited treatment of line mix- 
ing (aka line coupling) effects and from the assumed line pro- 
file. Line mixing in LBLRTM is modeled using a first order 
pertur bation approach. For CO2, line mixing is treated accord- 
ing to lNiro et all d2005h . LBLRTM employs a Voigt line shape 
at all atmospheric levels with an algorith m based on a linea r 
combination of approximating functions (IClough et all |2005). 
A Rautian, Galatry or speed-dependent Voigt profile would be 
needed to account more precisely for spee d-dependent broad- 
ening effects and line mixing effects (e.g., iBrault et"aTL [2003, 
Fig.l). The errors induced by this simplification are, however, 
comparatively small against the uncertainties in other line pa- 
rameters entries of the HITRAN database. 



4.2.3. Limitations from the model atmosphere 

As long as the total column densities of the different molecular 
species can be varied, the uncertainties in the model atmosphere 
have only a minor impact on the model spectrum. Exceptions 
are CO and O3, where changes in the vertical abundance pro- 
files lead to notably different line widths, because either the 
abundance maximum or minimum of the species is close to the 
temperature inversion of the Earths atmosphere (see Fig. Q] and 
The model for our mid-infrared case presented in Sec. 14.1.51 
does show a dependence on the type of model atmosphere, i.e. 
whether an equatorial or night-time MIPAS model is used. The 
intrinsic line width is, however, degenerate with the instrumental 
profile and while the formal resolution measured as the width of 
the instrumental profile varies, the residuals give no hint which 
of the two model atmospheres is more representative of the ac- 
tual conditions. 

For the water vapour, we did not find a significant improve- 
ment in the fit quality when using the MM5 atmospheric mod- 
els instead of the GDAS models. The fitting delivered the same 
PWV for both models and the different vertical distribution of 
the water vapour had no noticeable influence on our test cases. 
Re-fitting of the vertical distribution of the water vapour will 
thus lead only to minor improvements, if at all, and seems only 
feasible when enough lines are present, preferably from different 
vibrational transitions. 



4.2.4. Limitations from stellar lines 

The fitting of a telluric model to an observed spectrum will be 
more complicated in the presence of blends between telluric and 
stellar lines. Various strategies could be used to allow for the 
same fitting quality as compared to the featureless early type 
stars used in this study. (I) The usage of a stellar model spectrum 
in the fitting process would significantly improve the telluric 
model. (II) In the thermal infrared, the telluric emission lines 
could be used to further constrain the fit in terms of wavelength 
solution and molecular abundances. (Ill) In the worst case, an 
observation of a telluric standard star in the same spectral set- 
ting but not necessarily at the same time and air mass can be used 
to co nstrain the model fit for the target (e.g., IVillanueva et all 
l2008h . 

5. Conclusions 

We have shown that with the use of LBLRTM, a line-by-line ra- 
diative transfer code, and an appropriate model atmosphere, syn- 
thetic telluric transmission spectra can be constructed that match 
the observed telluric spectrum to 2% or better (see, e.g., Fig|5j, 
even close to saturated lines. These results were achieved with 
only a modest effort in back-fitting of atmospheric properties. 
Compilation of LBLRTM is straight forward and when using 
the described IDL fitting routines for x 2 minimisation, the time 
needed for a model fit to a single CRIRES chip is less than 2 min- 
utes on standard desktop PC . User interaction is only required 
to download the appropriate GDAS model from the NOAA web- 
site, to set up an initial wavelength scale, and to choose the 
species to be fitted. 

Although full photon noise limited performance is never 
achieved, the usage of the synthetic transmission spectra offers 
an alternative to the usage of standard star observations. The 
tight constraints in matching airmass and instrumental profile of 
the science target imposed on telluric standard stars are hard to 
meet in reality without downgrading the achievable SNR or the 
spectral resolution. 

We demonstrate this problem by using the data shown in 
Sect. 14.1.21 this time showing in FigQ~T]the spectrum from nod- 
ding position A divided by the one in nodding position B, taken 
1.5 min later. The difference in airmass is only 0.005. The spec- 
trum of the B position had to be mapped onto the wavelength 
scale of the A position due to the slit curvature which changes 
the wavelength zero point and the dispersion between both posi- 
tions. We note that at this step the wavelength solution provided 
by the telluric model for both nodding positions had to be used 
to achieve this result. Despite the fact that both spectra are of the 
same star and taken close in time and airmass, a slight change 
in the AO performance over time and in the instrumental profile 
between both nodding positions led to the apparent residuals in 
the division. The detector response at the two nodding positions 
is uniform to better than 0.5% after correcting for non-linearity 
and applying a flatfield correction as outlined in Sec. 2. Hence, 
changes in the instrumental profile are the reason for the ob- 
served residuals. 

Under more realistic circumstances, i.e. when a different tar- 
get is observed after a longer time interval and not matching the 
airmass that well, the expected residuals from a simple division 
by another telluric standard star would be even larger. Even in 
this very favorable case, the residuals in the line cores are at 
the order of 10% and the residuals in the wings are comparable 
to the ones after dividing through a telluric model (see Fig. |5j- 
Acquiring accurate empirical telluric spectra with CRIRES is 
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Table 2. Precipitable water vapour values (in mm PWV at zenith). 



Date & Time 


GDAS 




MM5 




ESO a 


CUT) 


predicted 


fitted 


predicted 


fitted 


predicted 


2006-10-10 03:00 


1.8 


0.8 


1.8 


0.8 


1.9 


2007-03-01 06:00 


5.3 


2.6 






3.1 


2007-03-03 09:00 


4.1 


1.6 


2.4 


1.6 


1.2 


2007-03-05 06:00 


10.8 


8.6 






3.9 


2007-04-22 09:00 


9.6 


7.2 






2.5 


2007-05-10 01:00 


2.0 


0.4" 






3.0 


2008-02-29 09:00 


3.6 


2.2 






2.8 



a Prediction for Paranal Observatory, based on satellite measurements. 

Obtained from http: //www. eso . org/gen-fac/pubs/astclim/forecast/meteo/ERASMUS/ 
b PWV determined from water vapour lines not shown in Fig. 1101 



challenging and the usage of telluric model spectra is an inter- 
esting alternative to the observation of standard stars, especially 
since it also saves valuable observing time. 

An effective removal of telluric absorption lines demands, 
but in most cases also provides, a wavelength calibration at an 
extremely high level of precision. Narrow and steep absorption 
lines are sensitive to wavelength errors at the level of 1/50 pixel 
(30 m/s). Having the telluric absorption lines directly imprinted 
in the science spectrum, allows for the most unbiased wave- 
length calibration when compared to telluric emission lines or 
calibration lamp lines since the telluric absorption lines experi- 
ence the same IP as the stellar lines. Given the large variation in 
accuracy of the line positions in the HITRAN database, a wave- 
length solution based on a weighted fit to telluric lines offers fur- 
ther improvements, especially when a sufficient number of lines 
is present. With the excellent line positions of the 12 CC>2 lines 
shown in Fig.[5]a wavelength solution with an r.m.s. of < 20 m/s 
can be achieved. 

To a certain extent, the Earth's atmosphere serves as 
a 'gascelF for precise wavelength calibration as shown by 
ISeifahrt & Kaufll (12008) who used N 2 lines at 4100nm and 
reached a short term wavelength precision of < 10 m/s over sev- 
eral hours. 



1 T T T 1 , , , 1 
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2058 2060 2062 2064 2066 2068 

Fig. 11. Division of two spectra from consecutive nodding posi- 
tions of the same star shown in Fig. [5] The spectrum from nod- 
ding position B was mapped on the wavelength scale of of the 
spectrum from nodding position A before the division. The dif- 
ference in airmass is only 0.005. See Fig. |5]and text for details. 



Recently Figu eira et alj (1201 Oh have demonstrated, based on 
6 years of HARPS optical data, that atmospheric O2 absorp- 
tion lines are stable down to 10 m/s and can be used as a fre- 
quency standard with an internal precision of better than 2 m/s. 
Transferring this precision into the infrared wavelength regime 
needs a precise model for the telluric interference spectrum to 
account for variability, mostly H2O content and airmass. The lat- 
ter is not trivial in the presence of fully or nearly fully saturated 
lines. 

The technique described in the current paper was used to 
calculate telluric model spectra as part of the data an alysis of a 
near In frared radial velocity program conducted by iBean et aTl 
d2010 a b) with the CRIRES spectrograph. Used in conjunction 
with an Ammonia gascell, an RV precision of 3-5 m/s was ob- 
tained on mid-M dwarfs in a spectral region near 2300 nm. 

Despite not being intended as a study for time dependent 
water vapour distribution, we find that the observed variability 
in the PWV above Paranal, ranging from 0.4 mm to nearly 9 mm 
in our observations, exceeds the range of satellite based predic- 
tions for the PWV by far. We have compiled all predicted and 
fitted PWV values used in this study in Tab. [2] The PWV can 
change notably in the course of hours and re-fitting of the to- 
tal PWV for a specific observation is unavoidable. A detailed 
study on the water vapour column density above Paranal will be 
presented elsewhere (Smette et al., in preparation). The GDAS 
meteorological model very often predicts too high PWV values, 
which we have always scaled globally, without changing its ver- 
tical distribution. Comparing GDAS and MM5 models for the 
same nights, results in the same PWV value and the use of the 
more realistic MM5 model did not significantly improve or fit 
results. 

Choosing well characterised lines, we do see the potential in 
the CRIRES data to retrieve atmospheric abundance information 
for atmospheric trace gases at or below the 1 % level delivering 
important information about the short and long term variations 
of greenhouse gases from astronomical observations. Such mea- 
surements could be a valuable byproduct of astronomical spectra 
collected with high-resolution infrared spectrographs. 
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